{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd4e8894-a336-4cc1-ae65-abd9e6ec0378",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f7fd41-500b-489f-a9a0-15a9508a0fba",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "55f0b0f1-91ef-427b-8980-af46954e588f",
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "cannot import name '_lazywhere' from 'scipy._lib._util' (/home/sc/miniconda3/envs/py3.13/lib/python3.13/site-packages/scipy/_lib/_util.py)",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mImportError\u001b[39m                               Traceback (most recent call last)",
      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[3]\u001b[39m\u001b[32m, line 6\u001b[39m\n\u001b[32m      4\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mseaborn\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01msns\u001b[39;00m\n\u001b[32m      5\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mscipy\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m stats\n\u001b[32m----> \u001b[39m\u001b[32m6\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mformula\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mapi\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m ols, glm\n\u001b[32m      7\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mstats\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01manova\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m anova_lm\n\u001b[32m      8\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mstats\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mmulticomp\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m pairwise_tukeyhsd\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/py3.13/lib/python3.13/site-packages/statsmodels/formula/api.py:2\u001b[39m\n\u001b[32m      1\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mregression\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mlinear_model\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mlm_\u001b[39;00m\n\u001b[32m----> \u001b[39m\u001b[32m2\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mdiscrete\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mdiscrete_model\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mdm_\u001b[39;00m\n\u001b[32m      3\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mdiscrete\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mconditional_models\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mdcm_\u001b[39;00m\n\u001b[32m      4\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mregression\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mmixed_linear_model\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mmlm_\u001b[39;00m\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/py3.13/lib/python3.13/site-packages/statsmodels/discrete/discrete_model.py:38\u001b[39m\n\u001b[32m     36\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mbase\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01m_parameter_inference\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mpinfer\u001b[39;00m\n\u001b[32m     37\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mbase\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m _prediction_inference \u001b[38;5;28;01mas\u001b[39;00m pred\n\u001b[32m---> \u001b[39m\u001b[32m38\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mdistributions\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m genpoisson_p\n\u001b[32m     39\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mregression\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mlinear_model\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mlm\u001b[39;00m\n\u001b[32m     40\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mtools\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m data \u001b[38;5;28;01mas\u001b[39;00m data_tools, tools\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/py3.13/lib/python3.13/site-packages/statsmodels/distributions/__init__.py:7\u001b[39m\n\u001b[32m      2\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01m.\u001b[39;00m\u001b[34;01mempirical_distribution\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m (\n\u001b[32m      3\u001b[39m     ECDF, ECDFDiscrete, monotone_fn_inverter, StepFunction\n\u001b[32m      4\u001b[39m     )\n\u001b[32m      5\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01m.\u001b[39;00m\u001b[34;01medgeworth\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m ExpandedNormal\n\u001b[32m----> \u001b[39m\u001b[32m7\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01m.\u001b[39;00m\u001b[34;01mdiscrete\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m (\n\u001b[32m      8\u001b[39m     genpoisson_p, zipoisson, zigenpoisson, zinegbin,\n\u001b[32m      9\u001b[39m     )\n\u001b[32m     11\u001b[39m __all__ = [\n\u001b[32m     12\u001b[39m     \u001b[33m'\u001b[39m\u001b[33mECDF\u001b[39m\u001b[33m'\u001b[39m,\n\u001b[32m     13\u001b[39m     \u001b[33m'\u001b[39m\u001b[33mECDFDiscrete\u001b[39m\u001b[33m'\u001b[39m,\n\u001b[32m   (...)\u001b[39m\u001b[32m     21\u001b[39m     \u001b[33m'\u001b[39m\u001b[33mzipoisson\u001b[39m\u001b[33m'\u001b[39m\n\u001b[32m     22\u001b[39m     ]\n\u001b[32m     24\u001b[39m test = PytestTester()\n",
      "\u001b[36mFile \u001b[39m\u001b[32m~/miniconda3/envs/py3.13/lib/python3.13/site-packages/statsmodels/distributions/discrete.py:5\u001b[39m\n\u001b[32m      3\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mscipy\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mstats\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m rv_discrete, poisson, nbinom\n\u001b[32m      4\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mscipy\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mspecial\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m gammaln\n\u001b[32m----> \u001b[39m\u001b[32m5\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mscipy\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01m_lib\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01m_util\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m _lazywhere\n\u001b[32m      7\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mstatsmodels\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mbase\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mmodel\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m GenericLikelihoodModel\n\u001b[32m     10\u001b[39m \u001b[38;5;28;01mclass\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mgenpoisson_p_gen\u001b[39;00m(rv_discrete):\n",
      "\u001b[31mImportError\u001b[39m: cannot import name '_lazywhere' from 'scipy._lib._util' (/home/sc/miniconda3/envs/py3.13/lib/python3.13/site-packages/scipy/_lib/_util.py)"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy import stats\n",
    "from statsmodels.formula.api import ols, glm\n",
    "from statsmodels.stats.anova import anova_lm\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n",
    "from statsmodels.discrete.count_model import ZeroInflatedNegativeBinomialP, ZeroInflatedPoisson\n",
    "from patsy import dmatrix\n",
    "# import statsmodels.api as sm\n",
    "from hmmlearn import hmm\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "\n",
    "# 尝试导入statsmodels，处理兼容性问题\n",
    "try:\n",
    "    import statsmodels.api as sm\n",
    "    use_statsmodels = True\n",
    "except ImportError as e:\n",
    "    if '_lazywhere' in str(e):\n",
    "        print(\"警告: statsmodels与当前SciPy版本不兼容，将跳过回归分析部分\")\n",
    "        use_statsmodels = False\n",
    "    else:\n",
    "        raise\n",
    "\n",
    "\n",
    "\n",
    "# 设置随机种子保证结果可复现\n",
    "np.random.seed(42)\n",
    "\n",
    "# ========================\n",
    "# 1. 模拟车险数据集\n",
    "# ========================\n",
    "print(\"生成模拟车险数据...\")\n",
    "n = 5000  # 保单数量\n",
    "\n",
    "# 生成基本保单信息\n",
    "data = pd.DataFrame({\n",
    "    'policy_id': range(1, n+1),\n",
    "    'age': np.random.randint(18, 70, n),  # 车主年龄\n",
    "    'vehicle_value': np.random.lognormal(mean=3.0, sigma=0.5, size=n),  # 车辆价值\n",
    "    'driving_experience': np.random.randint(1, 50, n),  # 驾龄\n",
    "    'vehicle_type': np.random.choice(['Sedan', 'SUV', 'Truck', 'Sports'], n, p=[0.4, 0.3, 0.2, 0.1]),  # 车型\n",
    "    'vehicle_color': np.random.choice(['Red', 'Blue', 'Black', 'White', 'Silver'], n),  # 颜色\n",
    "    'area': np.random.choice(['Urban', 'Suburban', 'Rural'], n, p=[0.5, 0.3, 0.2]),  # 地区\n",
    "    'annual_mileage': np.random.gamma(shape=2, scale=5000, size=n),  # 年行驶里程\n",
    "    'credit_score': np.random.normal(loc=700, scale=50, size=n)  # 信用评分\n",
    "})\n",
    "\n",
    "# 添加索赔次数（使用零膨胀负二项分布）\n",
    "# 红色车辆索赔频率高8.2%，乡村地区索赔频率低12.3%\n",
    "red_effect = np.where(data['vehicle_color'] == 'Red', 1.082, 1.0)\n",
    "rural_effect = np.where(data['area'] == 'Rural', 0.877, 1.0)\n",
    "mean_claim = 0.15 * red_effect * rural_effect * (data['annual_mileage']/10000)\n",
    "claim_counts = np.random.negative_binomial(n=2, p=1/(1+mean_claim), size=n)\n",
    "zero_mask = np.random.random(n) < 0.65  # 65%的保单无索赔\n",
    "claim_counts[zero_mask] = 0\n",
    "data['claim_count'] = claim_counts\n",
    "\n",
    "# 添加索赔金额（使用伽玛分布）\n",
    "claim_amounts = np.random.gamma(shape=1.5, scale=2000, size=n)\n",
    "claim_amounts[claim_counts == 0] = 0  # 无索赔时金额为0\n",
    "data['claim_amount'] = claim_amounts\n",
    "\n",
    "# 添加高风险驾驶行为标记\n",
    "data['high_risk'] = np.where((data['annual_mileage'] > 20000) & \n",
    "                             (data['age'] < 25) & \n",
    "                             (data['vehicle_type'] == 'Sports'), 1, 0)\n",
    "\n",
    "print(\"数据样例：\")\n",
    "print(data.head())\n",
    "print(\"\\n数据描述：\")\n",
    "print(data.describe())\n",
    "\n",
    "# ========================\n",
    "# 2. 抽样分布分析\n",
    "# ========================\n",
    "print(\"\\n分析索赔次数的分布特征...\")\n",
    "\n",
    "# 可视化索赔次数分布\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(1, 2, 1)\n",
    "sns.histplot(data['claim_count'], kde=False, bins=15)\n",
    "plt.title('索赔次数分布')\n",
    "plt.xlabel('索赔次数')\n",
    "\n",
    "# 零膨胀特征分析\n",
    "zero_claim_ratio = (data['claim_count'] == 0).mean()\n",
    "print(f\"零索赔保单比例: {zero_claim_ratio:.2%}\")\n",
    "\n",
    "# 拟合零膨胀负二项分布 (ZINB)\n",
    "print(\"\\n拟合零膨胀负二项分布(ZINB)...\")\n",
    "exog = data[['age', 'vehicle_value', 'driving_experience', 'annual_mileage', 'credit_score']]\n",
    "exog = sm.add_constant(exog)  # 添加常数项\n",
    "\n",
    "# 创建并拟合ZINB模型\n",
    "zinb_model = ZeroInflatedNegativeBinomialP(\n",
    "    endog=data['claim_count'],\n",
    "    exog=exog,\n",
    "    exog_infl=exog,\n",
    "    inflation='logit'\n",
    ")\n",
    "zinb_results = zinb_model.fit(maxiter=100)\n",
    "print(zinb_results.summary())\n",
    "\n",
    "# ========================\n",
    "# 3. 假设检验\n",
    "# ========================\n",
    "print(\"\\n进行假设检验...\")\n",
    "\n",
    "# 检验1: 红色车辆是否索赔频率更高\n",
    "red_claims = data[data['vehicle_color'] == 'Red']['claim_count']\n",
    "other_claims = data[data['vehicle_color'] != 'Red']['claim_count']\n",
    "\n",
    "# 独立样本t检验\n",
    "t_stat, p_value = stats.ttest_ind(red_claims, other_claims, equal_var=False)\n",
    "print(f\"红色车辆索赔频率检验: t = {t_stat:.3f}, p = {p_value:.4f}\")\n",
    "if p_value < 0.05:\n",
    "    red_mean = red_claims.mean()\n",
    "    other_mean = other_claims.mean()\n",
    "    print(f\"红色车辆平均索赔次数: {red_mean:.3f}, 其他颜色: {other_mean:.3f}\")\n",
    "    print(f\"红色车辆索赔频率高 {(red_mean/other_mean-1)*100:.1f}%\")\n",
    "\n",
    "# 检验2: 地区间索赔频率差异 (ANOVA)\n",
    "print(\"\\n地区间索赔频率差异分析:\")\n",
    "area_anova = ols('claim_count ~ area', data=data).fit()\n",
    "anova_table = anova_lm(area_anova)\n",
    "print(anova_table)\n",
    "\n",
    "# 事后检验 (Tukey HSD)\n",
    "tukey_results = pairwise_tukeyhsd(data['claim_count'], data['area'], alpha=0.05)\n",
    "print(tukey_results)\n",
    "\n",
    "# ========================\n",
    "# 4. 多因素分析\n",
    "# ========================\n",
    "print(\"\\n进行多因素分析...\")\n",
    "\n",
    "# 4.1 索赔频率预测 (GAM模型)\n",
    "print(\"构建广义可加模型(GAM)预测索赔频率...\")\n",
    "\n",
    "# 创建样条项\n",
    "data['mileage_spline'] = dmatrix(\"bs(annual_mileage, df=4, degree=3)\", \n",
    "                                {\"annual_mileage\": data['annual_mileage']}, \n",
    "                                return_type='dataframe')\n",
    "\n",
    "# 构建GAM模型公式\n",
    "gam_formula = \"claim_count ~ age + driving_experience + C(vehicle_type) + \" \\\n",
    "              \"C(vehicle_color) + C(area) + mileage_spline\"\n",
    "\n",
    "# 拟合负二项回归模型 (用于计数数据)\n",
    "gam_model = glm(gam_formula, data=data, \n",
    "               family=sm.families.NegativeBinomial()).fit()\n",
    "print(gam_model.summary())\n",
    "\n",
    "# 4.2 关键因素分析\n",
    "print(\"\\n关键因素影响分析:\")\n",
    "\n",
    "# 计算变量重要性\n",
    "coefs = gam_model.params\n",
    "std_errors = gam_model.bse\n",
    "z_scores = coefs / std_errors\n",
    "p_values = gam_model.pvalues\n",
    "\n",
    "# 创建结果DataFrame\n",
    "factor_importance = pd.DataFrame({\n",
    "    'Variable': coefs.index,\n",
    "    'Coefficient': coefs.values,\n",
    "    'Std_Error': std_errors.values,\n",
    "    'z_score': z_scores.values,\n",
    "    'p_value': p_values.values\n",
    "})\n",
    "\n",
    "# 过滤显著变量\n",
    "significant_factors = factor_importance[\n",
    "    (factor_importance['p_value'] < 0.05) & \n",
    "    (~factor_importance['Variable'].str.contains('Intercept'))\n",
    "].sort_values('z_score', key=abs, ascending=False)\n",
    "\n",
    "print(\"显著影响因素:\")\n",
    "print(significant_factors[['Variable', 'Coefficient', 'p_value']].head(10))\n",
    "\n",
    "# 可视化关键因素影响\n",
    "plt.figure(figsize=(10, 6))\n",
    "sns.barplot(x='Coefficient', y='Variable', \n",
    "           data=significant_factors.head(10), palette='viridis')\n",
    "plt.title('车险索赔频率的关键影响因素')\n",
    "plt.xlabel('回归系数')\n",
    "plt.ylabel('变量')\n",
    "plt.tight_layout()\n",
    "\n",
    "# ========================\n",
    "# 5. 随机过程 (驾驶行为建模)\n",
    "# ========================\n",
    "print(\"\\n使用隐马尔可夫模型(HMM)建模驾驶行为...\")\n",
    "\n",
    "# 模拟驾驶行为数据\n",
    "n_days = 365  # 一年数据\n",
    "n_drivers = 100\n",
    "\n",
    "# 驾驶状态: 0=安全, 1=中等风险, 2=高风险\n",
    "states = ['Safe', 'Medium', 'High']\n",
    "\n",
    "# 转移概率矩阵\n",
    "transition_matrix = np.array([\n",
    "    [0.85, 0.12, 0.03],  # 安全状态下的转移概率\n",
    "    [0.20, 0.70, 0.10],  # 中等风险\n",
    "    [0.05, 0.15, 0.80]   # 高风险\n",
    "])\n",
    "\n",
    "# 观测指标: 急刹车次数、超速事件、夜间驾驶\n",
    "emission_means = np.array([\n",
    "    [0.2, 0.1, 0.3],  # 安全状态\n",
    "    [1.0, 0.8, 1.5],  # 中等风险\n",
    "    [3.0, 2.5, 4.0]   # 高风险\n",
    "])\n",
    "\n",
    "# 生成驾驶行为序列\n",
    "driver_data = []\n",
    "for driver in range(n_drivers):\n",
    "    # 随机初始状态\n",
    "    current_state = np.random.choice(3)\n",
    "    daily_behavior = []\n",
    "    \n",
    "    for day in range(n_days):\n",
    "        # 记录当前状态\n",
    "        daily_behavior.append({\n",
    "            'driver_id': driver,\n",
    "            'day': day,\n",
    "            'state': states[current_state],\n",
    "            'hard_braking': np.random.poisson(emission_means[current_state, 0]),\n",
    "            'speeding': np.random.poisson(emission_means[current_state, 1]),\n",
    "            'night_driving': np.random.poisson(emission_means[current_state, 2])\n",
    "        })\n",
    "        \n",
    "        # 转移到下一天状态\n",
    "        current_state = np.random.choice(3, p=transition_matrix[current_state])\n",
    "    \n",
    "    driver_data.extend(daily_behavior)\n",
    "\n",
    "driver_df = pd.DataFrame(driver_data)\n",
    "\n",
    "# 训练HMM模型\n",
    "print(\"训练隐马尔可夫模型...\")\n",
    "observations = driver_df[['hard_braking', 'speeding', 'night_driving']].values\n",
    "\n",
    "# 标准化观测值\n",
    "scaler = StandardScaler()\n",
    "observations_scaled = scaler.fit_transform(observations)\n",
    "\n",
    "# 创建并训练HMM\n",
    "model = hmm.GaussianHMM(n_components=3, covariance_type=\"diag\", n_iter=100)\n",
    "model.fit(observations_scaled)\n",
    "\n",
    "# 解码隐藏状态\n",
    "hidden_states = model.predict(observations_scaled)\n",
    "driver_df['predicted_state'] = [states[int(s)] for s in hidden_states]\n",
    "\n",
    "# 评估模型效果\n",
    "accuracy = (driver_df['state'] == driver_df['predicted_state']).mean()\n",
    "print(f\"状态预测准确率: {accuracy:.2%}\")\n",
    "\n",
    "# 可视化状态转移\n",
    "plt.figure(figsize=(10, 6))\n",
    "state_counts = driver_df.groupby('day')['predicted_state'].value_counts(normalize=True).unstack()\n",
    "state_counts.plot.area(stacked=True, colormap='viridis')\n",
    "plt.title('驾驶行为状态随时间变化')\n",
    "plt.ylabel('状态比例')\n",
    "plt.xlabel('天数')\n",
    "plt.legend(title='驾驶状态')\n",
    "plt.tight_layout()\n",
    "\n",
    "# ========================\n",
    "# 6. 综合应用：动态保费定价\n",
    "# ========================\n",
    "print(\"\\n构建综合定价模型...\")\n",
    "\n",
    "# 合并索赔数据和驾驶行为\n",
    "driver_claim_data = data.sample(n_drivers).reset_index(drop=True)\n",
    "driver_claim_data['driver_id'] = range(n_drivers)\n",
    "merged_data = pd.merge(driver_df, driver_claim_data, on='driver_id')\n",
    "\n",
    "# 构建定价模型\n",
    "pricing_model = glm(\n",
    "    \"claim_count ~ predicted_state + age + vehicle_value + annual_mileage\",\n",
    "    data=merged_data,\n",
    "    family=sm.families.NegativeBinomial()\n",
    ").fit()\n",
    "\n",
    "# 预测索赔频率\n",
    "merged_data['predicted_claim'] = pricing_model.predict()\n",
    "\n",
    "# 计算基础保费\n",
    "base_premium = 1000  # 基础保费\n",
    "merged_data['premium'] = base_premium * (1 + merged_data['predicted_claim'] * 0.5)\n",
    "\n",
    "# 高风险驾驶附加费\n",
    "risk_surcharge = {\n",
    "    'Safe': 0.0,\n",
    "    'Medium': 0.15,\n",
    "    'High': 0.40\n",
    "}\n",
    "merged_data['risk_surcharge'] = merged_data['predicted_state'].map(risk_surcharge)\n",
    "merged_data['total_premium'] = merged_data['premium'] * (1 + merged_data['risk_surcharge'])\n",
    "\n",
    "# 分析不同状态下的保费差异\n",
    "premium_by_state = merged_data.groupby('predicted_state')['total_premium'].mean()\n",
    "print(\"\\n不同驾驶状态的平均保费:\")\n",
    "print(premium_by_state)\n",
    "\n",
    "# 保存结果\n",
    "merged_data.to_csv('auto_insurance_pricing.csv', index=False)\n",
    "print(\"分析完成! 结果已保存到 auto_insurance_pricing.csv\")\n",
    "\n",
    "# 显示所有图表\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a639fed-203e-4389-984f-c28c645a1b15",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "518d8dce-661a-4567-89b3-55d324492fd3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1a24eb7-2bfb-4474-9cd5-4fed0401dc86",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30ff84ef-4421-447d-96b5-f42c0cc4d7c2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "878c2795-d3fb-4c05-b403-fb1846a1e887",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8921640-f9bc-4140-b2df-2943e64f9907",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
